#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time    : 2021/3/28 18:24
# @Author  : Rem~
# @File    : 005SineSignal_WithNoise.py
# @function: 在004的基础上，添加噪声

import numpy as np
import matplotlib.pyplot as plt

# sample frequency
fs = 100
# signal frequency
f = 5
# signal phase
fai = 55*np.pi/180
# sample point
N = 1000
# x axis
x = np.linspace(0, N-1, N)
# amplitude
amp = 2**12
bias = amp
# 探测器接收的光信号
sig = amp*np.cos(2*np.pi*f*(x/fs)+fai)
# 实际带偏置的光信号
sig_bias = sig+bias
# 高斯白噪声，需要调用随机数
noise_amp = 2**4
noise = np.random.uniform(-noise_amp, noise_amp, N)

sig_bias_noise = sig_bias+noise
# 高斯窗函数
wn = np.exp(-((x-N/2)/N*4)**2)
# 高斯窗下含噪声的模拟光信号
sig_bias_noise_wn = sig_bias_noise*wn


fig = plt.figure("Time&Frequency Domain Analyse")
axes1 = fig.add_subplot(3, 2, 1)
axes2 = fig.add_subplot(3, 2, 2)
axes3 = fig.add_subplot(3, 1, 2)
axes4 = fig.add_subplot(3, 2, 5)
axes6 = fig.add_subplot(3, 2, 6)

axes1.plot(x, sig_bias_noise, "-b")
axes1.set(xlabel="sample point", ylabel="Time Domainsignal amplitude(LSB)")
axes1.grid(True)

axes2.plot(wn)
axes2.set(xlabel="sample point", ylabel="Gauss Window amplitude(LSB)")
axes2.grid(True)

axes3.plot(sig_bias_noise_wn)
axes3.set(xlabel="sample point", ylabel="Signal in Gauss Window amplitude(LSB)")
axes3.grid(True)

sig_bias_noise_wn_fft = np.fft.fft(sig_bias_noise_wn)
sig_bias_noise_wn_fft_magnitude = np.abs(sig_bias_noise_wn_fft)/N*2
sig_bias_noise_wn_fft_phase = np.angle(sig_bias_noise_wn_fft)*(180/np.pi)

axes4.stem(x[0:int(N/2)-1]*(fs/N), sig_bias_noise_wn_fft_magnitude[0:int(N/2)-1])
axes4.set(xlabel="frequency axis(MHz)", ylabel="signal amplitude(LSB)")
axes4.grid(True)

axes6.plot(x[0:int(N/2)-1]*(fs/N), sig_bias_noise_wn_fft_phase[0:int(N/2)-1])
axes6.set(xlabel="frequency axis(MHz)", ylabel="signal phase(°)")
axes6.grid(True)


sig_bias_noise_wn_fft_abs_max = np.max(sig_bias_noise_wn_fft_magnitude[5:int(N/2)-1])
sig_bias_noise_wn_fft_abs_max_f = np.argmax(sig_bias_noise_wn_fft_magnitude[5:int(N/2)-1])+5
sig_bias_noise_wn_fft_max_f_phase = sig_bias_noise_wn_fft_phase[sig_bias_noise_wn_fft_abs_max_f]

print('频谱峰值:', sig_bias_noise_wn_fft_abs_max,
      '\n最大峰值点:', sig_bias_noise_wn_fft_abs_max_f*(fs/N), 'MHz',
      '\n最大频率点相位', sig_bias_noise_wn_fft_max_f_phase, '°')

plt.show()
